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Abstract  identification.  Fogel  [2J  has  proposed  the  use  ol'  EP 

in  system  identification  for  evolving  ARM  AX  pr<x;e.s.s 
Evolutionary  programming  (EP)  has  been  coefficients.  His  investigations  extend  to 

demonstrated  to  be  an  effective  method  of  system  simultaneous  parameter  and  model  structure 

identification  of  single-input-single-output  (SISO)  estimation.  'rhi.s  paper  extends  the  work  on  .SiSO 

systems.  This  paper  investigates  the  use  of  EP  in  systems  to  SIMO  systems.  The  extension  to  multi¬ 
system  identification  of  single-input-multi -output  input,  multi-output  (MIMO)  .systems  is  obvious. 

(SIMO)  systems.  EP  is  used  to  identify  parameters  Identification  of  a  SIMO  system  could  be 

of  a  linear,  time-invariant  system.  Specifically,  this  divided  into  identifying  a  set  of  SISO  system  transfer 

paper  examines  the  identification  of  SIMO  systems  functions  between  the  input  and  each  of  the  outputs, 

whose  measurements  caucoti'am  different  amounts  of  However,  each  of  these  transfer  functions  would  be 

noise.  A  cost  function  is  projxised  to  take  into  identified  using  only  the  measurements  of  its  output 

account  disparate  noisy  observations.  and  not  the  full  set  of  oEservations  which  contains 

more  information.  Also,  it  is  not  clear  how  one 
I.  Introduction  might  combine  these  individual  tran.sfer  functions 

with  parameter  estimates  of  varying  degrees  of 
System  identification  is  the  process  of  uncertainty  into  one  combined  system  model  that 

developing  an  accurate  model  of  a  dynamic  system.  would  be  required  to  control  ail  the  states  of  the 

This  involves  selecting  the  appropriate  model  system. 

structure  and  choosing  the  model  parameters  which  The  remainder  of  this  paper  is  organized  as 

optimize  an  objective  function.  The  usual  approach  follows.  Section  II  contains  an  overview  of  the  EP 

taken  in  system  identification  is  to  assume  a  model  paradigm.  Section  HI  formulates  EP  for  SIMO 

structure  based  on  the  process  physics  and  then  to  sy.stem  identification  and  de.scribes  the  model 

determine  the  model  parameters.  Oftentimes,  structure.  Section  IV  presents  a  linear,  time-invariant 

nonlinear  models  are  linearized  so  that  an  appropriate  system  that  is  u.sed  for  rompnter  expenments. 

state  space  formulation  may  be  used.  Section  V  di.scusses  the  results  and  Section  VI 

Typically,  methods  of  system  identification  provides  a  summary  of  the  results  and  conclusions, 

involve  minimizing  the  objective  function  using 

gradient  descent  techniques.  Use  of  the  gradient  may  El.  Evolutionary  Programming 

result  in  the  identified  parameters  being  only  locally 

optimal.  Continuity  constraints  are  another  limitation  EP  is  a  parallel  stocha.stic  search  technique 

associated  with  pure  gradient  search  techniques.  A  that  maps  natural  evolution  to  an  effective  multi-agent 

Kalman  filter  approach  is  given  in  [IJ  along  with  search  strategy.  This  .search  technique  is  not  ba.sed 

cautions  on  its  application  to  system  identification.  on  gradient  methods  and  is  therefore  not  su.sceptible 

Many  algorithms  have  been  proposed  which  attempt  to  entrapment  in  local  minima.  Search  by  simulated 

to  alleviate  the  inherent  problems  associated  with  evolution  wa-s  first  described  by  Eogel  (.Ij.  Ihis 


gradient  methods.  Recent  approaches  have  included 
simulated  annealing,  genetic  algorithms,  neural 
networks  and  evolutionary  programming.  The  latter 
method  ha.s  proven  extremely  successful  for  function 
optimization. 

Evolutionary  programming  (EP)  has  been 
succe.s.sfully  applied  to  the  problem  of  .system 


investigation  follows  the  .search  strategy  oiitiine  l  m  <j 
Fogel  [2]  and  described  by  the  sequence  below 
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Create  the  initial  population  P  consisting  of 
N  parents 

As.sc-ss  the  fitness  ot  each  parent  p,  in  the 
ponulalion. 
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Figure  1.  Evolutionary  programming 
parameter  estimation  block  diagram. 


3.  Mutate  each  element  in  the  population  by  a 
N(0,a)  random  variable,  where  a  is 
proportional  to  the  fitness. 

4.  Assign  a  new  fitness  score  to  each  element 
of  the  population. 

5.  Hold  a  stochastic  competiticn  between  the 
members  of  the  population. 

6.  Rank  the  members  of  the  population 
according  to  the  results  (number  of  "wins") 
of  step  5. 

7.  Repeat  step  3  with  the  highest  ranked  N 
members  from  the  population. 

This  algorithm  is  applied  to  system  identification  in 
the  following  section. 

ni.  EP  for  System  Identification 

The  linear  time-invariant  SIMO  structure  is 
assumed  to  be  of  the  form 

a:(A:+I)  =  Fx{k)+Gu{k) 
y(k)  =  Hx(k)+v{k) 

where  x(k)  is  the  state  vector  of  length  n,  u(k)  is  the 
input,  y(k)  is  the  vector  of  measurements  of  length  1, 
v(k)  is  the  vector  of  ii.t.dourei.icnl  noi.se.  F  is 

the  .system  matrix,  G  i.s  the  input  weighting  matrix 


ami  H  is  the  measurement  scaling  matrix.  The  bkxk 
diagram  for  estimating  parameters  asing  EP  is  shown 
in  Figure  1. 

The  vectors  to  be  evolved  will  consist  of  the 
elements  of  the  F  and  G  matnce.s.  The  H 
(measurement)  matrix  will  not  be  evolved  because 
finding  all  three  of  the  matrices  is  an  underdefined 
problem  since  an  infinite  number  of  combinations  of 
F.  G  and  H  matrices  apply  to  the  same  system 
dynamics  For  this  rea.son,  the  measurements  will  be 
assigned  to  the  first  /  states  of  the  system.  This 
results  in  an  //  matrix  of  the  form 

1 

Given  this  formulation,  the  method  of 
applying  EP  to  the  identification  problem  is 
straightforward.  The  remaining  issues  to  be  resolved 
are  the  mutation  strategy  and  the  objective  function. 

Two  mutation  strategies  were  investigated. 
The  first  incorporates  a  standard  deviation  of  the 
mutation  for  each  vector  that  is  linearly  proportional 
to  the  fitness  score  as  suggested  by  Fogel  [2].  This 
results  in  the  size  of  the  mutations  decreasing  as  the 
population  approaches  the  optimum  value.  In 
addition  to  ensuring  convergence,  this  strategy  also 
assures  an  effective  search  by  allowing  larger 
mutations  for  elements  that  have  poor  fitne.ss  scores. 
Atmar  [4]  suggests  a  variety  of  mechanism.s  to 
accelerate  the  search  process. 

^  The  second  mutation  strategy  inve.stigated  is 
based'on  the  number  of  iterations.  The  mutation  of 
any  element  in  the  vector  can  be  described  by  the 
perturbation 

5  =  Me'""* 

where  is  a  N(0,a)  random  variable  (r.v.),  r  is  a 
U[0,1I  r.v.,  a  is  a  constant  and  k  is  the  generation 
number.  This  mutation  is  a  function  of  the  number 
of  iterations  only.  A  random  walk  is  maintained  over 
a  high  number  of  generations  while  still  preserving 
the  capability  of  infinitesimal  perturbations  to  achieve 
optimum  values  as  k  increa,se.s.  Reasonable  results 
were  obtained  using  this  mutation  strategy  tor 
optimizing  mobile  manipulator  configurations  (.‘i). 

The  objective  function  must  be  carefully 
chosen  for  this  search.  Not  only  will  the 
measurements  generally  be  of  different  units,  hui  ilk- 
amount  ot  confidence  in  each  measurement  may 
differ  due  to  sensor  characteristics.  .A  cost  lunction 


has  been  formulated  to  take  into  account  disparate 
sensor  noise  and  the  difference  between  state 
estimates  and  observations.  The  error  vector  between 
observation  and  estimate  is  given  by 

y  =  y-y 


where  y  is  the  observation  and  y  is  the  estimate.  The 
cost  function  incorporates  sensor  uncertainty  (noise) 
using  the  diagonal  standard  deviation  matrix 


Q  = 


0  0 


0 


0  0  . 


a 


where  a,  corresponds  to  the  noise  characteristics  of 
the  sensor.  The  cost  associated  with  each 
parameter  estimate  can  be  described  as 

r  = 

n 


IV.  Simulated  System 

A  simple  second-order  system  based  on 
spnng-ma.ss-damper  dynamic.s  was  u.sed  for  computer 
expenments.  Parameters  were  chosen  so  that  the 
sy.stem  was  cntically  damped.  The  continuous  state 
et}uations  can  be  written  as 


r 

j 

‘  - 

0  1 

1 

0 

X  - 

-1  -2 

,t  +■  : 

1 

y  = 


l  0 
0  1 


X  *■  y 


The  equivalent  discretized  state-space  formulation, 
assuming  a  zero-order-hold  and  a  sample  rate  of 
T=0.05  sec,  is  given  by 


0.9988  0.0476 
-0.0476  0.9037 


x(k)  + 


0.0012 

0.0476 


Kk) 


y{k)  = 


I  0 
0  1 


x(.k)  +  i>(k) 


where  the  variances  are  selected  according  to  sensor 
specifications,  n  is  the  length  of  the  error  vector  and 
R  ■■=  n^O  or  the  identity  matrix  I  depending  on  the 
experiment.  Using  R  =  accomplishes  two 
things.  First,  it  allows  the  cost  to  be  a  dimensionless 
quantity.  Second,  it  scales  each  measurement  by  the 
amount  of  confidence  in  its  accuracy.  If  the 
measurement  is  noisy,  its  contribution  to  the  cost 
function  is  reduced  since  the  prediction  errors  could 
be  due  primarily  to  measurement  noise.  If  an 
observation  is  relatively  noise-free,  its  contribution  to 
the  cost  function  will  be  large  since  the  error  is  due 
primarily  to  estimation  errors,  assuming  the  model 
structure  is  correct.  For  these  trials,  it  i.s  assumed 
that  the  variances  are  constant  but  not  necessarily 
equal.  If  the  variances  are  equal,  then  the  cost 
function  is  simply  scaled  by  this  magnitude.  Using 
R  =  I  and  scaling  the  cost  function  will  yield  the 
same  ie.sults. 


The  forcing  function  is  ensured  to  be  persistently 
exciting  by  defining  u(k)  as  an  r.v.  with  a  uniform 
distribution.  The  system  was  propagated  forward 
100-200  time  steps  for  each  simulation  from  an  initial 
state  randomly  chosen  from  a  Ul-1,1]  distribution. 

V.  Results 

Simulations  were  conducted  using  both 
mutation  strategies.  In  general,  the  second  mutation 
strategy  outperformed  the  first.  The  reason  for  this 
is  .still  under  investigation.  The  state  variables  had  to 
be  constrained  to  prevent  numeric  overflow.  This  is 
not  unreasonable  given  that  the  time  integration  of  an 
unstable  sy.stem  is  not  bounded.  However,  it  doe.s 
have  implications  in  system  identification  of  unstable 
plants. 

The  first  expenmenl  tor  the  second  mutation 
strategy  u,sed  numerically  equal  measurement  nf^i-se 
lenns  wUh  N’(0, 0.0.5)  siati.stic.s  and  if-0.01.  The 
re,sulting  parameters  are  given  in  Table  I  under 
experiment  I.  The  average  percentage  error  for  all 


parameters  is  approximately  8%. 

The  magnitude  of  the  forcing  function  was 
increased  tenfold  in  an  attempt  to  obtain  better  system 
identification.  The  parameters  resulting  from  this 
experiment  are  given  in  Table  1  in  the  experiment  2 
column.  As  expected,  these  parameters  are  have  a 
lower  overall  average  error  than  those  obtained  in 
experiment  1.  For  this  case,  the  average  overall 
parameter  error  is  roughly  3.14%. 

Figure  2  shows  the  mean  and  best  population 
scores  using  the  second  mutation  strategy.  The  cost 
function  is  simply  the  mean  sum  squared  error  (MSB) 
of  both  observations.  At  roughly  300  iterations  the 
identification  process  converged.  This  run  had 
equivalent  noise  statistics  on  each  measurement. 
Figure  3  illustrates  the  convergence  of  the  first 
parameter  (fj,)  in  the  highest  ranked  F  matrix.  The 
other  parameters  converge  in  a  similar  fashion. 

An  experiment  was  conducted  to  compare 
the  effect  of  the  two  cost  functions  using  the  second 
mutation  strategy.  The  standard  deviation  of  the 
position  measurement  noise  was  0.01  and  the 
standard  deviation  of  the  velocity  measurement  noise 
was  0.05.  fn  both  cases  a  was  set  at  0.05.  The  best 
values  of  the  parameters  are  shown  in  Table  I.  The 
parameter  values  under  the  experiment  3  heading  are 
for  the  cost  function  that  included  the  variance  term 
{R  =  O’^Q).  The  parameter  values  given  under  the 
experiment  4  heading  are  for  R  equal  to  the  identity 
matrix.  The  experiment  that  incorporated  a  cost 
function  without  the  variance  term  unexpectedly  had 
better  results  than  the  experiment  which  incorporated 
a  confidence  weighting  term,  The  reasons  for  this 
result  are  still  being  investigated. 

VI.  Conclusions 

EP  is  a  powerful  search  technique  that  can 
be  applied  to  sy.stem  identification.  This  research  has 
demonstrated  that  EP  can  be  extended  to  identify 
linear  SIMO  systems  of  known  order.  The  extension 
to  identification  of  MIMO  systems  is  obvious.  It  is 
al.so  clear  that  the  scaling  of  the  cost  function  must  be 
chosen  carefully  to  account  for  disparate  sen.sor  data 
and  measurement  noi.se.  The  results  given  in  this 
work  have  yet  to  be  statistically  verified  and  therefore 
should  not  be  generalized.  Investigations  on  the 
algorithmic  properties  of  EP  are  being  continued  for 
validation  in  further  experiments. 

The  formulation  de.scribed  here  is  a  general 
one  for  the  identificahuii  of  state-.space  models,  but 
there  are  still  additional  considerations  to  be 


addres,sed.  First,  the  mode!  ignored  process  noi.se, 
the  terms  that  account  for  differences  from  the  actual 
dynamic  equations.  Second,  this  work  assumed  that 
the  order  of  the  sy.stem  is  known.  To  identify  the 
order  of  the  model  as  well  as  the  parameters,  a  multi¬ 
dimensional  equivalent  of  Akaike’.s  information 
critenon  (AIC)  or  minimum  defining  length  (MDL) 
pnnciple  similar  to  that  used  by  Fogel  [6]  for  neural 
networks  must  be  developed. 

The  application  of  EP  for  adaptive  control  is 
readily  apparent.  Issues  to  be  addres.sed  in  using  BP 
for  adaptive  control  pertain  to  the  bandwidth  of  the 
process  being  controlled.  That  is,  of  course, 
assuming  a  traditional  adaptive  control  architecture. 
EP  could  be  implemented  in  an  arrangement  such  that 
system  identification  would  not  have  to  be  explicit. 
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Figure  2.  The  mean  and  minimum  (best)  Figure  3.  Evolution  of  parameter  f„. 

population  costs. 


Para¬ 

meters 

Actual 

values 

Final  Best  Value  Statistics 

1**  Experiment 

2“*  Experiment 

3"*  Experiment 

4’*'  Experiment 

value 

%  error 

value 

%  error 

value 

%erTor 

value 

%enor 

f.. 

0.9988 

0.9965 

0.23 

0.9986 

0.02 

1.0193 

2.05 

0.9973 

0.15 

f,2 

0.0476 

0.0527 

10.7 

0.0480 

0.84 

0.1015 

113.2 

0.0422 

11.3 

I  — 

-0.0476 

-0.0449 

5.67 

-0.0481 

HQH 

-0.0745 

56.51 

-0.504 

5.88 

^22 

0.9031 

0.90163 

0.23 

0.9033 

0.04 

0.8315 

7.99 

0.9062 

0.28 

gi 

0.0012 

0.0009 

25.0 

0.0010 

16.7 

-0.0152 

1367 

0.0009 

25.0 

g2 

0.0476 

0.0507 

6.51 

0.0477 

0.21 

0.0480 

0.84 

0.0440 

7.56 

Table  I.  Comparison  of  final  value  results  of  experiments  and  actual  parameters. 


